Stretching skin: The physiological limit and beyond

aDepartment of Mechanical Engineering, Stanford University, Stanford, CA 94305, USA

bDepartment of Plastic Surgery, Rainbow Babies and Children’s Hospital, Case Western Reserve University, Cleveland, OH 44106, USA

cDepartments of Mechanical Engineering, Bioengineering, and Cardiothoracic Surgery, Stanford University, Stanford, CA 94305, USA

The goal of this manuscript is to establish a novel computational model for skin to characterize its constitutive behavior when stretched within and beyond its physiological limits. Within the physiological regime, skin displays a reversible, highly nonlinear, stretch locking, and anisotropic behavior. We model these characteristics using a transversely isotropic chain network model composed of eight wormlike chains. Beyond the physiological limit, skin undergoes an irreversible area growth triggered through mechanical stretch. We model skin growth as a transversely isotropic process characterized through a single internal variable, the scalar-valued growth multiplier. To discretize the evolution of growth in time, we apply an unconditionally stable, implicit Euler backward scheme. To discretize it in space, we utilize the finite element method. For maximum algorithmic efficiency and optimal convergence, we suggest an inner Newton iteration to locally update the growth multiplier at each integration point. This iteration is embedded within an outer Newton iteration to globally update the deformation at each finite element node. To illustrate the characteristic features of skin growth, we first compare the two simple model problems of displacement- and force-driven growth. Then, we model the process of stretch-induced skin growth during tissue expansion. In particular, we compare the spatio-temporal evolution of stress, strain, and area gain for four commonly available tissue expander geometries. We believe that the proposed model has the potential to open new avenues in reconstructive surgery and rationalize critical process parameters in tissue expansion, such as expander geometry, expander size, expander placement, and inflation timing.

Skin is the largest organ of our body. It covers a surface area of about two square meters and weighs up to four kilograms. Skin is a composite material consisting of three layers, a thin protective outer layer, the epidermis, a thick elastic inner layer, the dermis, and a subcutaneous base layer of fatty tissue, the hypodermis [60]. Made up of a complex network of collagen and elastin fibers, it is mainly the dermis which gives the skin its mechanical strength [11]. Within its physiological limits, skin almost behaves like rubber: Its mechanical response is highly nonlinear [51], initially weak, but much stiffer at higher stretch levels, limited through a characteristic locking stretch [52]. In contrast to rubber though, skin is highly anisotropic [37], with a larger stiffness along pronounced collagen fiber orientations which manifests itself macroscopically in the form of Langer’s lines [38]. When chronically stretched beyond its physiological limit, skin displays a fascinating behavior: It increases its surface area to reduce the mechanical load [18].

The controlled surface area growth through mechanical stretch was first proposed more than half a century ago to reconstruct a traumatic ear [50]. Tissue expansion has since then revolutionized reconstructive skin surgery and is now widely used to repair birth defects [7, 25], correct burn injuries [4, 17], and reconstruct breasts after tumor removal [54]. Tissue expansion is the perfect strategy to grow skin that matches the color, texture, and hair bearance of the surrounding healthy skin, while minimizing scars and risk of rejection [56]. Figure 1 illustrates an example of tissue expansion in pediatric forehead reconstruction [26]. The patient, a one-year old boy, presented with a giant congenital nevus concerning 25 percent of the forehead, affecting the hairline and the left cheek. Approximately one in 20,000 infants is born with giant congenital nevi, typically creating significant aesthetic distortion to the involved anatomical site, associated with severe psychological distress to the patients and their families [6]. In the patient shown in Figure 1, simultaneous forehead, cheek, and scalp expanders were used for in situ tissue growth. Tissue expander placement is typically performed via dissection of a subcutaneous pockets adjacent to the skin defect, while the expander ports to regulate expander filling are either buried in a remote location, or left outside for ease of injection. The amount of filling is controlled by visual inspection of skin color, capillary refill, and simple palpation of the stretched skin [56]. Multiple subsequent serial inflations stretch the skin and stimulate tissue growth. Once new skin is produced, the device is removed, and the new skin is used to repair the adjacent defect zone. The follow-up photograph in Figure 1, right, shows the patient at age three after forehead reconstruction. Similar expansion techniques have successfully been used to grow skin in the trunk [5], and in the upper and lower extremities [25].

Tissue expansion for pediatric forehead reconstruction. The patient, a one-year old boy presented with a giant congenital nevus concerning approximately 25 percent of the forehead, affecting the hairline and the cheek. Simultaneous forehead, cheek, and scalp expanders were implanted for in situ skin growth. This technique allows to resurface large anatomical areas with skin of similar color, quality, and texture. The follow-up photograph shows the patient at age three after forehead reconstruction.

Figure 2 shows a schematic sequence of the mechanical processes that occur during tissue expansion. Initially, at biological equilibrium, the skin is in a natural state of resting tension [60], left. Tissue expanders are implanted subcutaneously and gradually inflated, stretching the skin beyond the physiological limit, top. This triggers a series of stretch-induced signaling pathways [65]. Mechanotransduction affects a network of integrated cascades including cellular architecture and function such as cytoskeletal structure, extracellular matrix, enzyme activity, second messenger systems, and ion channel activity [18]. As a consequence, the skin grows and restores the state of resting tension, right.

Schematic sequence of tissue expander inflation. At biological equilibrium, the skin is in a physiological state of resting tension, left. A tissue expander is implanted subcutaneously between the skin, consisting of the epidermis and dermis, and the hypodermis. When the expander is inflated, the skin is loaded in tension, top. Mechanical stretch induces cell proliferation causing the skin to grow. Growth restores the state of resting tension, right.

This cycle of expander inflation, overstretch, growth, and relaxation is repeated multiple times, typically on a weekly basis, throughout a period of 6–8 weeks [16, 41]. Remarkably, as demonstrated by immunocytochemical analyses, the expanded tissue undergoes normal cell differentiation and maintains its characteristic phenotype [71]. Although the tissue initially displays epidermal thickening and dermal thinning upon expansion, both thickness changes are usually reversible under expander removal [68].

Several studies have focused on understanding the adaptation of skin from a biomechanical point of view. A process parameter that has received major attention is the geometry of the tissue expander [14]. Figure 3 displays four commonly used tissue expander geometries, circular, square, rectangular, and crescent-shaped. For regular circular and rectangular expanders, simple mathematical models have been proposed to kinematically correlate the expander volume to its surface area [19, 59]. From an engineering point, it is intuitive though that purely kinematic models severely overestimate the net gain in surface area [69]. With a discrepancy of up to a factor four, these models assume that the entire deformation can be attributed to irreversible growth, completely neglecting the elastic deformation which is reversible upon expander removal [41]. In an attempt to account for this error, empirical correction factors of 6.00, 3.75, and 4.50 have been proposed for circular, rectangular, and crescent-shaped expanders [69]. Despite these efforts, the choice of the appropriate expander geometry and size is still almost exclusively based on the surgeon’s personal preference, and the discrepancy between recommended shapes, sizes, and volumes remains enormous [41]. This demonstrates the ongoing need to rationalize criteria for a standardized device selection.

Tissue expanders to grow skin flaps for defect correction in reconstructive surgery. Typical applications are birth defect correction, scar revision in burn injuries, and breast reconstruction after tumor removal. Devices are available in different shapes and sizes, circular, square, rectangular, and crescent-shaped. They consist of a silicone elastomer inflatable expander with a reinforced base for directional expansion, and a remote silicone elastomer injection dome. Reprinted with permission, Mentor Worldwide LLC.

In this manuscript, we propose a rigorous, mechanistic approach to systematically compare different tissue expander geometries in terms of stress, strain, and area gain. To model skin growth in response to tissue expansion, we adopt the framework of finite growth based on the multiplicative decomposition of the deformation gradient into an elastic and a growth part [57]. A tremendous amount of research has been devoted to establish continuum theories for finite growth within the last decade [20, 42]. These theories have been applied successfully to characterize growing tumors [2], tendons [21], vascular tissue [35, 64], cardiac tissue [22, 55]. While earlier studies were primarily of theoretical and analytical nature [10, 63], we can now observe a clear trend towards the computational modeling of volumetric growth, typically by introducing the growth tensor as an internal variable within a finite element framework [27, 58]. We would like to point out, however, that these theories, although successful in characterizing growth on a macroscopic tissue level, remain phenomenological in nature. We would have to consult more sophisticated mixture theories [1, 29], if we wanted to understand the microstructural origin of growth. A recent monograph that compares different approaches to growth summarizes the essential findings, trends, and open questions in this progressively evolving new field [3].

Despite ongoing research in growing biological systems, the growth of thin biological membranes remains severely understudied [24]. Motivated by a first study on axisymmetric skin growth [62], we have recently proposed a prototype model for growing isotropic membranes to model skin expansion in a general three-dimensional setting [15]. The goal of this manuscript is to significantly refine this initial isotropic model and incorporate the basic features of skin including its extreme nonlinearity, its locking behavior, and its anisotropic nature [33, 34], to precisely quantify stress, strain, and area gain in response to different, arbitrarily shaped tissue expander geometries.

This manuscript is organized as follows. In Section 2, we give a brief overview of the continuum equations for finite growth including the kinematic equations, the balance equations, and the constitutive equations. In Section 3, we illustrate the temporal and spatial discretizations of the biological and mechanical equilibrium equations, along with their consistent algorithmic linearizations. We then demonstrate the features of our model in Section 4 by comparing the simple model problems of displacement- and force-driven growth, and illustrate skin growth in response to tissue expander inflation. We close with a brief discussion, some limitations, and concluding remarks in Section 5.

In this Section, we illustrate the governing equations for skin growth which consist of three basic sets of equations; in Section 2.1 the kinematic equations based on the concept of an incompatible growth configuration and the multiplicative decomposition of the deformation gradient; in Section 2.2 the balance equations of open systems phrased in their mass specific format; and in Sections 2.3 and 2.4 the constitutive equations, both for overstretch-induced growth and for the baseline elastic response.

We adopt the kinematics of finite deformations and introduce the deformation map ϕ, which, at any given time t ∈ , maps the material placement X of a physical particle in the material configuration to its spatial placement x in the spatial configuration .

In what follows, we apply a formulation which is entirely related to the material frame of reference. Accordingly, ∇{○} = ∂X {○}|t and Div {○} = ∂X {○}|t: G−1 denote the gradient and the divergence of any field {○}(X,t) with respect to the material placement X at fixed time t. Herein, G−1 is the contravariant material metric. To characterize finite growth, we introduce an incompatible growth configuration [40, 61], and adopt the multiplicative decomposition of the linear tangent map [57],

into a reversible elastic part Fe: T → T and an irreverisble growth part Fg: T → T , see Figure 4. This implies that the total Jacobian

Kinematics of finite growth. Illustration of covariant spatial metric g, deformation tensors C and Ce, stress tensors S, P, Pe and τ, and mappings F = Fe · Fg and F−t = Fe −t · Fg −t between tangent spaces T and cotangent spaces T * in the material configuration, the intermediate configuration, and the spatial configuration [27, 48].

obeys a similar multiplicative decomposition into an elastic part Je = det (Fe) and a growth part Jg = det (Fg). We idealize skin as a thin layer characterized through the unit normal n0 in the undeformed reference configuration. The length of the deformed skin normal n = cof(F) · n0 = J F−t · n0 introduces the area stretch,

which we can again decompose into an elastic area stretch ϑe = ||cof(Fe) · ng/||ng|| || and a growth area stretch ϑg = ||cof(Fg) · n0 ||, where ng = cof(Fg) · n0 = Jg Fg − t · n0 denotes the grown skin normal. Here, cof(○) = det(○) (○)−t denotes the cofactor of the second order tensor (○), such that ϑ = J [n0 · C−1 · n0]1/2. As characteristic deformation measures, we introduce the right Cauchy Green tensor C and its elastic counterpart Ce as the pull backs of the spatial metric g to the undeformed reference configuration and to the intermediate configuration,

where both are related through the following identity, Ce = Fg − t · C · Fg − 1. In the context of finite growth, we typically parameterize the constitutive equations in terms of the reversible elastic deformation tensor Ce. Its relevant invariants

and their derivatives

then take the above representations, where we have introduced the direction of transverse isotropy ν0 characterizing the collagen orientation in the skin layer, or, macroscopically speaking, through the direction of Langer’s lines [38]. Finally, we introduce the pull back of the spatial velocity gradient l to the intermediate configuration,

which obeys the additive split into the elastic velocity gradient Le = Fe−1 ·Ḟe and the growth velocity gradient Lg = Ḟg ·Fg−1. Here, we have applied the notation {○̇} = ∂t{○}|X to denote the material time derivative of any field {○}(X,t) at fixed material placement X. Figure 4 illustrates the kinematics of finite growth in terms of the covariant spatial metric g, the deformation tensors C and Ce, and the mappings F = Fe · Fg and F−t = Fe − t · Fg − t between tangent and cotangent spaces T and T * in the material configuration, the intermediate configuration, and the spatial configuration.

We characterize the growing skin layer using the framework of open system thermodynamics in which the material density ρ0 is allowed to change as a consequence of growth [30, 32]. The balance of mass for open systems balances its rate of change ρ̇0 with a possible in- or outflux of mass R and mass source [53, 67].

Similarly, the balance of linear momentum balances the density-weighted rate of change of the momentum v̇, where v = ϕ̇ is nothing but the spatial velocity, with the momentum flux P = F · S and the momentum source ρ0 b,

here stated in its mass-specific form [31, 32]. P and S are the first and second Piola-Kirchhoff stress tensors, respectively. Last, we would like to point out that the dissipation inequality for open systems

typically contains an extra entropy source ρ0 to account for the growing nature of living biological systems [30, 44]. Equations (10) and (11) represent the mass-specific versions of the balance of momentum and of the dissipation inequality [31].

To characterize the growing tissue, we introduce the constitutive equations for the mass flux R, for the mass source , and for the growth tensor Fg. In what follows, we assume that mass changes induced by diffusion are significantly smaller than local changes in mass. Accordingly, we assume that the mass flux R is negligible,

and that all changes in mass can be attributed exclusively to the mass source . Immunocytochemistry has shown that expanded tissue undergoes normal cell differentiation [71]. Accordingly, we assume that the newly grown skin has the same density and microstructure as the initial tissue. This implies that the mass source

can be expressed as the density-weighted trace of the growth velocity gradient tr (Lg) = Ḟg: Fg − t [27]. Motivated by clinical observations [56], we represent growth as a strain-driven, transversely isotropic, irreversible process, characterized through a single growth multiplier ϑg that reflects the irreversible area increase perpendicular to the skin normal n0.

For this particular transversely isotropic definition of the growth tensor, for which the material is not allowed to grow in the thickness direction [68], the area growth is identical to the volume growth, i.e., ϑg = det(Fg) = Jg. Because of the simple rank-one update structure, we can apply the Sherman-Morrison formula to invert the growth tensor in explicit form.

It introduces the following simple expression for the growth velocity gradient, Lg=ϑ.g/ϑgI+[1−ϑ.g/ϑg]n0⊗n0 which proves convenient to explicitly evaluate the mass source as R0=ρ0[1+2ϑ.g/ϑg]. Motivated by physiological observations of stretch induced skin expansion, we introduce a strain-driven evolution law for the growth multiplier.

To control unbounded growth, we introduce the weighting function

where τ denotes the adaptation speed, γ calibrates the shape of the adaptation curve, and ϑmax denotes the maximum area growth [27, 42]. The growth criterion

is driven by the elastic area stretch ϑe = ϑ/ϑg, such that growth is activated only if the elastic area stretch exceeds a critical physiological stretch limit ϑcrit, where 〈○〉 denote the Macaulay brackets.

Once we have evaluated the irreversible part of the model, i.e., the growth multiplier ϑg, the growth tensor Fg, the growth velocity gradient Lg, and the mass source , we can turn to evaluate the reversible elastic part of the model, i.e., the momentum flux S and the momentum source b. For the sake of simplicity, we assume the latter to vanish identically, b = 0. We model skin as a transversely isotropic elastic material that can be characterized through the Helmholtz free energy ψ = ψ̂ (C, Fg, ν0 ⊗ ν0), parameterized in terms of the right Cauchy Green tensor C, the growth tensor Fg, and the preferred material orientation ν0 characterizing the direction of Langer’s lines. We adopt a transversely isotropic eight chain model for which the free energy of a representative eight chain unit cell consists of three contributions [34] as illustrated in Figure 5.

Transversely isotropic eight chain model. Individual chains are modeled as wormlike chains with an energy ψwlc parameterized in terms of the end-to-end length r. Eight chains are assembled in a transversely isotropic unit cell with dimensions a and b, and a characteristic orientation ν0. The energy of each unit cell consists of the bulk energy ψblk, the energy of the eight individual chains ψchn, and their repulsive contributions ψrep.

The first term ψblk is purely isotropic and captures the effect of bulk incompressibility in terms of the third invariant I3e [21]. The second term ψchn reflects the effective assembly of the eight individual chain energies and introduces a constitutive coupling between the first and fourth invariants I1e and I4e [33]. The third term ψrep is the repulsive term that accounts for an initial stress-free reference configuration [12]. For the individual chains, we adopt a wormlike chain model based on the single chain energy,

where A is the persistence length, L is the contour length, r is the end-to-end length of the chain, θ is the absolute temperature, and k is the Boltzmann constant [13]. In a transversely isotropic unit cell with dimensions a and b, the undeformed end-to-end length is r0=a2+2b2/2, and the deformed end-to-end length r is a function of the first and fourth invariant I1 and I4 as introduced in equation (6).

Based on these considerations, we can introduce the individual energy terms,

where we have used the following abbreviation for the repulsive weighting factor ψ̄rep.

Here, γblk and γchn are the chain and bulk densities, and β is a macroscopic bulk parameter. Using the free energy (19), we can now evaluate the dissipation inequality (11).

Similar to finite strain plasticity, we observe that the Mandel stress of the intermediate configuration Me = Ce · Se is energetically conjugate to the growth velocity gradient Lg = Ḟg · Fg − 1. Moreover, from the dissipation inequality (24), we obtain the definition of the second Piola Kirchhoff stress S as thermodynamically conjugate quantity to the right Cauchy Green deformation tensor C.

Here, we have introduced the second Piola Kirchhoff stress of the intermediate configuration

in terms of the individual stress contributions corresponding to the three energy terms introduced in equations (22).

Here, we have introduced the following abbreviations for the second order bases of the chain stress and of the repulsive stress.

The basis of the chain stress is a result of the derivative of the end-to-end length r with respect to the elastic right Cauchy Green tensor, S̄chn = 8 r dr/dCe. The basis of the repulsive stress S̄rep = 2d ψ̄rep/dCe is used to construct the repulsive energy ψ̄rep in equation (23) such that the initial state r = r0 is stress free, Srep(r0) ≐ −Schn(r0) and thus S̄rep(r0) ≐ S̄ chn(r0). Figure 6 displays the transversely isotropic nature of our eight chain model. The dots represent experimental measurements from uniaxial tests on rabbit skin tested parallel and perpendicular to Langer’s lines [12, 39]. The lines represent the corresponding computational simulation. The model nicely captures the characteristic features of skin, including the strong non-linearity, the anisotropy, and the locking stretches. Last, we can evaluate the elastic constitutive moduli Le by taking the second derivative of the Helmholtz free energy ψ with respect to the elastic right Cauchy Green tensor Ce.

Uniaxial tension test. Transversely isotropic wormlike-chain based eight chain model. Dots represent experimental measurements on rabbit skin tested parallel and perpendicular to Langer’s lines [39]. Lines represent the corresponding computational simulation. Doted lines represent the parallel and perpendicular locking stretches. The model nicely captures the characteristic features of skin, including the strong non-linearity, the anisotropy, and the locking stretches.

Its three individual contributions take the following forms.

The fourth order bases of the chain term L̄chn = S̄chn ⊗ S̄chn and of the repulsive term L̄rep = 2 dS̄rep/dCe can be expressed as follows.

In equation (30), denotes the fourth order material identity defined as = [G−1 ⊗̄ G−1 + G−1 ⊗ G−1]/2. Here, we have applied the abbreviations ⊗̄ and ⊗ for the non–standard dyadic products according to the following component-wise definitions {• ⊗̄○}i jkl = {•}ik ⊗ {○}jl and {•⊗○}i jkl = {•}il ⊗ {○}jk. Last, we would like to reiterate that the choice of the particular free energy (19) only indirectly affects the growth process itself. The growth model is inherently modular and can easily be combined with more or less complex baseline elasticity models [15].

The classical eight chain model [8] naturally follows as a special case of our general chain network model by assuming that the cell dimensions a and b take equal values.

This implies that the chain energy ψchn and the repulsive energy ψrep become functions of the first invariant I1e alone. The overall free energy ψ no longer depends on the fiber direction ν0, rather, it characterizes an isotropic chain network response [34, 49].

The special case transverse isotropy follows by assuming a degenerated unit cell for which b tends to zero.

For this model, all chains are oriented in a single direction ν0. The resulting fiber terms ψchn and ψrep are thus no longer functions of the first invariant I1e. They depend exclusively on the stretch of the chains represented through the fourth invariant I4e. However, this type of model, which has been termed decoupled reinforcement model [46, 47], completely neglects network effects which are a characteristic feature of skin.

The governing equations for finite growth introduced in the previous section are complex and highly nonlinear. In this section, we illustrate their computational solution within an incremental iterative nonlinear finite element framework. To characterize the growth process at each instant in time, we introduce the growth multiplier ϑg as an internal variable, and solve its evolution equation (16) locally at each integration point using the finite difference method. To explore the interplay between growth and mechanics, we discretize the governing equations for finite growth (2), (10), and (25) in space using the finite element method. In this section, we first derive the discrete local residual and the corresponding tangent moduli for the local Newton iteration to iteratively determine the growth multiplier ϑg. Then, we derive the stresses S for the discrete global residual and the constitutive moduli L for the iteration matrix of the global Newton iteration to iteratively determine the deformation ϕ.

To discretize the biological equilibrium equation (16) in time, we partition the time interval of interest into nstp subintervals,

and focus on the interval [tn, tn+1] for which Δt = tn+1 − tn > 0 denotes the current time increment. Our goal is to determine the current growth multiplier ϑg for a given deformation state F at time t, and a given growth multiplier ϑng at the end of the previous time step tn. For the sake of compactness, here and from now on, we omit the index (○)n+1 for all quantities at the end of the current time step tn+1. To approximate the material time derivative of the growth multiplier ϑ̇, we introduce the following finite difference approximation.

In the spirit of implicit time stepping schemes, we now reformulate the evolution equation (16) with the help of this finite difference approximation, introducing the discrete residual Rϑ in terms of the unknown growth multiplier ϑϑ.

We suggest to solve this nonlinear residual equation for the unknown growth multiplier using a local Newton iteration. Within each iteration step, we calculate the linearization of the residual Rϑ with respect to the growth multiplier ϑg.

From equations (17) and (18), we can extract the linearizations of the weighting function ∂kg/∂ϑg = −γkg/[ϑmax − ϑg] and of the growth criterion ∂φg/∂ϑg = −ϑ/ϑg2. Within each iteration step, we calculate the iterative update of the unknown growth multiplier ϑg ← ϑg − Rϑ/Kϑ until convergence is achieved, i.e., until the local growth update Δϑg = −Rϑ/Kϑ is below a user-defined threshold value. In what follows, we will assume negligible mass diffusion, R = 0. This implies that, if necessary, the remaining balance of mass, ρ.0=ρ0[1+2ϑ.g/ϑg] can simply be evaluated locally in a post-processing step once local convergence is achieved.

With the simplifying assumptions of a vanishing momentum source, b = 0, and negligible inertia effects, v̇ = 0, the mechanical equilibrium equation (10) reduces to the internal force balance, Div (F · S) ≐ 0. We cast it into its weak form, ∇δϕ: [F · S] dV ≐ 0, through the multiplication with the test function δϕ and the integration over the domain of interest , to solve it globally on the node point level. To discretize it in space, we partition the domain of interest into nel finite elements B0e.

Our goal is to determine the deformation state ϕ for a given loading at time t. To approximate the test function d δϕ, the unknown deformation ϕ, and their gradients ∇δϕ and ∇ϕ, we apply an isoparametric Bubnov-Galerkin based finite element interpolation,

where Ni,N j are the element shape functions and i, j = 1,..,nen are the element nodes. We now reformulate the weak form of the balance of linear momentum (10) with the help of these finite element approximations, introducing the discrete residual RIϕ in terms of the unknown nodal deformation ϕJ.

Herein, the operator A symbolizes the assembly of all element residuals at the j = 1,..,nen element nodes to the global residual at the global node points J = 1,..,nel. We can evaluate the global discrete residual (38), once we have iteratively determined the growth multiplier ϑg for the given deformation state F and the given history ϑng as described in Section 3.1. Then we can successively determine the growth tensor Fg from equation (14), the elastic tensor Fe = F · Fg − 1 from equation (2), the elastic stress Se from equation (25), and lastly, the second Piola Kirchhoff stress S in terms of the stress Se in the intermediate configuration (26).

Again, we suggest an incremental iterative Newton algorithm to solve the nonlinear residual equation for the unknown deformation (38). The consistent linearization of the residual RIϕ with respect to the nodal vector of unknowns ϕJ introduces the global stiffness matrix.

The fourth order tensor L denotes the Lagrangian constitutive moduli which we can determine directly from the linearization of the Piola Kirchhoff stress S with respect to the total right Cauchy Green tensor C.

The first term

represents nothing but the pull back of the elastic moduli Le onto the reference configuration, where Le = 2∂Se/∂Ce are the constitutive moduli of the elastic material model as introduced in equation (29). The second term

consists of two contributions that resemble a geometric and a material stiffness contribution known from nonlinear continuum mechanics. The third term

and the fourth term

depend on the particular choice for the growth tensor Fg in equation (14) and on the evolution equation for the growth multiplier ϑg in equation (16), respectively. For each global Newton iteration step, we iteratively update the current deformation state ϕ←ϕ−KIJϕ−1·RIϕ until we achieve algorithmic convergence. Upon convergence, we store the corresponding growth multipliers ϑg at the integration point level. To solve these nonlinear finite element equations, we implement the growth model in a custom-designed version of the multipurpose nonlinear finite element program FEAP [66].

In contrast to the uniaxial tension test in Figure 6 which is only showing the acute, purely reversible elastic response of stretched skin tissue, we now explore how skin would respond chronically if we stretched it beyond its physiological limits. First, we illustrate the conceptual features of stretch-induced growth by studying two simple model problems, displacement-and force-driven skin expansion of a simple prototype sheet. Then, we focus on the physiological problem of tissue expansion by exploring different expander geometries and their different placements with respect to Langer’s lines [38]. If not stated otherwise, we choose the elastic material parameters for skin, calibrated by means of the uniaxial tension experiments shown in Figure 6 [12, 39]. Accordingly, our contour length is A = 1.85, our persistence length is L = 2.125, our eight chain unit cell dimensions are a = 2.43 and b = 1.95 [33], our chain and bulk densities are γchn = 1.75 · 1021 and γblk = 100, and our macroscopic bulk parameter is β = 4.5 [21]. The Boltzmann constant is k = 1.30 · 10−23J/K, and the absolute temperature is θ = 310K. To characterize the growth process, we choose the maximum area growth to ϑmax = 2.4 [15], the elastic stretch limit to ϑcrit = 1.80 for the model problem in Section 4.1 and to ϑcrit = 1.12 for the tissue expansion problem in Section 4.2, the adaptation speed to τ = 1.0, and the shape parameter for the adaptation curve to γ = 2.0 [15]. For the sake of illustration, in the following examples, we normalize the expansion time to t = 1.0, noting that the growth process scales linearly in time with the adaptation speed τ. Common tissue expansion procedures involve repeated cycles of expander inflation, overstretch, growth, and relaxation, typically repeated on a weekly basis, throughout a period of multiple weeks [16, 41]. To clearly illustrate the influence of anisotropy in the following examples, we slightly decrease the cell dimension to b = 2/3 a in comparison to the calibrated response of Figure 6 where b = 4/5 a.

Let us first illustrate the conceptual features of our stretch-induced growth model by exploring the two simple model problems of displacement- and force-driven skin expansion of a square skin sheet of size 1.0 ×1.0 and thickness 0.2.

For the displacement-driven case, we linearly increase the prescribed displacements from 0.0 at time t = 0.0 to 1.0 at time t = 0.4. This implies that the skin sheet is stretched biaxially to a final size of 2.0 × 2.0, i.e., total area stretch is ϑ = 4.0. The displacements are then held constant from t = 0.4 to t = 1.0 to allow the tissue to grow. Figure 7 illustrates the resulting temporal evolution of the total area stretch ϑ, the reversible elastic area stretch ϑe, and the irreversible area growth ϑg. The horizontal dashed lines represent the elastic stretch limit beyond which skin growth is activated ϑcrit, and the maximum area growth ϑmax. Displacements are increased linearly up to the vertical dashed line at t = 0.4. At this time, the total area stretch is ϑ = 4.000, the elastic area stretch ϑe = 2.157, and the area growth ϑg = 1.855. Then, displacements are held constant to allow the skin to grow until t = 1.0. At this time, growth has almost completely converged to the biological equilibrium state with a total area stretch of ϑ = 4.000, an elastic area stretch of ϑe = 1.866, and an area growth of ϑg = 2.143. The curves confirm, that, at all times, the multiplicative decomposition of the deformation gradient F = Fe ·Fg introduced in equation (2) carries over to the multiplicative decomposition of the total area stretch ϑ = ϑe ϑe of equation (4). Obviously, upon displacement control, we observe material relaxation indicated through a gradual increase of growth at a constant total stretch, while the elastic stretch, and, accordingly the stresses, decrease.

Relaxation test. Temporal evolution of total area stretch ϑ, reversible elastic area stretch ϑe, and irreversible growth area stretch ϑg for displacement driven skin expansion. Horizontal dashed lines represent the elastic stretch limit beyond which skin growth is activated ϑcrit, and the maximum area growth ϑmax. Displacements are increased linearly up to the vertical dashed line and then held constant. Displacement control induces relaxation indicated through the gradual decrease in elastic stretch and stress, while the growth stretch increases at a constant total stretch.

For the force-driven case, we linearly increase the prescribed forces from 0.0 at time t = 0.0 to 128.0 at time t = 0.4. The forces are then held constant from t = 0.4 to t = 1.0 to allow the tissue to grow. Figure 8 illustrates the corresponding temporal evolution of the total area stretch ϑ, the reversible elastic area stretch ϑe, and the irreversible growth area stretch ϑg. Again, the horizontal dashed lines represent the elastic stretch limit ϑcrit beyond which skin growth is activated, and the maximum area growth ϑmax. Forces are increased linearly up to the vertical dashed line at t = 0.4. At this time, the total area stretch is ϑ = 3.646, the elastic area stretch is ϑe = 2.048, and the area growth is ϑg = 1.780. Then, the forces are held constant to allow the skin to grow until t = 1.0. At this time, growth has almost completely converged to the biological equilibrium state with a total area stretch of ϑ = 4.557, an elastic area stretch of ϑe = 2.055, and an area growth of ϑg = 2.217. Obviously, upon force control, we observe material creep indicated through a gradual increase of growth and of the total stretch, while the elastic stretch and the stresses remain constant.

Creep test. Temporal evolution of total area stretch ϑ, reversible elastic area stretch ϑe, and irreversible growth area stretch ϑg for force driven skin expansion. Horizontal dashed lines represent the elastic stretch limit beyond which skin growth is activated ϑcrit, and the maximum area growth ϑmax. Forces are increased linearly up to the vertical dashed line and then held constant. Force control induces creep indicated through the gradual increase in growth stretch and total stretch at constant elastic stretch and stress.

Finally, we focus on the physiological problem of tissue expansion and compare different expander geometries and their different placements with respect to Langer’s lines. As illustrated in Figure 9, we model skin as a 0.2 cm thin 12 ×12 cm2 square sheet, discretized with 3 ×24 ×24 = 1728 trilinear brick elements, with 4 × 25 × 25 = 2500 nodes and 7500 degrees of freedom. To explore the impact of different tissue expander geometries, we model a circular, a square, a rectangular, and a crescent shaped expander, and place them in alignment with and orthogonal to Langer’s lines. For the sake of comparison, the base surface area of all four expanders is scaled to a size of A0 = 37 cm2. In tissue expansion, expander placement is performed via dissection of a subcutaneous pocket between the approximately 0.2 cm thick skin layer consisting of the epidermis and the dermis [60], and the thick fatty hypodermis. In Figure 9, the subcutaneous pocket is displayed in light red, while the intact, non-dissected tissue is displayed in white. As illustrated in Figure 9, we model the expansion procedure by pressuring the 148 light red elements of the subcutaneous pocket in which the expander is placed, while fixing the bottom nodes of all remaining white elements. We increase the pressure linearly up to a maximum pressure of p = 160 at t = 0.08. We then keep the pressure constant and watch the skin grow until we achieve convergence towards the biological equilibrium state at t = 1.00.

Tissue expansion and skin growth. Skin is modeled as a 0.2 cm thin 12 × 12 cm2 square sheet, discretized with 3 × 24 × 24 = 1728 trilinear brick elements, with 4 × 25 × 25 = 2500 nodes and 7500 degrees of freedom. We explore different tissue expanders, circular, square, rectangular, and crescent shaped, when placed in alignment with and orthogonal to Langer’s lines. The base surface area of all expanders is scaled to 148 elements corresponding to 37 cm2. Expanders are placed in a subcutaneous pocket between the skin and the hypodermis, here shown in light red, while the intact, non-dissected tissue is displayed in white. The light red area is gradually pressurized from underneath, while the bottom nodes of all intact white skin elements are fixed.

Figure 10 displays the temporal evolution of the fractional area gain for all four expander geometries. The circular expander displays the largest fractional area gain, followed by the rectangular expander aligned with Langer’s lines, the square expander, the crescent shaped expander aligned with Langer’s lines, and, lastly, the crescent shaped and rectangular expanders placed orthogonal to Langer’s lines. The curves demonstrate the characteristic creep-type growth under constant pressure, a response which is conceptually similar to the one discussed in Figure 8. The fractional area gain converges towards a biological equilibrium state after a characteristic time period, here characterized through the normalized time t = 1.0, in reality a period of several days. Tissue expansion is a gradual procedure that is performed through repeated cycles of expander inflation, overstretch, growth, and relaxation. Tissue growth stagnates after several days, and expander re-inflation becomes necessary. This is why tissue expansion patients need to see their physician on a weekly basis for tissue expander re-inflation. The overall procedure takes place throughout a period of multiple weeks [41].

Tissue expansion and skin growth. Temporal evolution of fractional area gain. The expander pressure is increased gradually up to the first vertical dashed line and then held constant to allow the skin to grow. The rectangular and crescent shaped expanders display a larger fractional area gain when aligned with Langer’s lines than when placed orthogonal to Langer’s lines. The curves demonstrate the characteristic creep-type growth under constant pressure with a gradual convergence towards the biological equilibrium state.

Figure 11 illustrates the temporal evolution of the corresponding expander volume. A constant pressure causes the tissue to creep and the expander volume to increase. Again, the curves clearly mimic the anisotropic response of skin. Under the same pressure, applied to the same area, the final expander volumes are larger for the rectangular and crescent-shaped expanders when aligned with Langer’s lines than when placed orthogonal to Langer’s lines.

Tissue expansion and skin growth. Temporal evolution of expander volume. The expander pressure is increased gradually up to the first vertical dashed line and then held constant to allow the skin to grow. The rectangular and crescent shaped expanders display a larger expander volumes when aligned with Langer’s lines than when placed orthogonal to Langer’s lines. The curves demonstrate the characteristic creep-type growth under constant pressure with a gradual convergence towards the biological equilibrium state.

Figure 12 demonstrates the spatio-temporal evolution of the area growth ϑg for the circular, square, rectangular, and crescent-shaped expanders. Each expander is placed in two different positions, in alignment with and orthogonal to Langer’s lines, as indicated through the vectors ν0. The corresponding growth contours are shown in the upper and lower row of each set. The color code ranges from ϑg = 1.0 for the initially ungrown state, shown in blue, to ϑg = ϑmax = 2.4 for the fully grown state, shown in red. Snapshots correspond to t = 0.08, t = 0.16, t = 0.24, t = 0.48, and t = 0.96, from left to right, corresponding to the vertical dashed lines in Figures 10 and ​and1111.

Tissue expansion and skin growth. Spatio-temporal evolution of growth area stretch ϑg for circular, square, rectangular, and crescent-shaped expanders with two different orientations with respect to Langer’s lines ν0 for each set. Growth is larger when the expanders are positioned along material’s strong direction, aligned with Langer’s lines. Under the same pressure applied to the same base surface area, the circular expander in rows 1 and 2 induces the largest amount of growth; the rectangular and crescent-shaped expanders placed orthogonal to Langer’s lines in rows 6 and 8 induce the smallest amount of growth. The color code illustrates the evolution of the growth multiplier ϑg, ranging from ϑg = 1.0 for the initially ungrown skin, shown in blue, to ϑg = ϑmax = 2.4 for the fully grown state, shown in red. Snapshots correspond to t = 0.08, t = 0.16, t = 0.24, t = 0.48, and t = 0.96, from left to right, corresponding to the labels in Figures 10 and ​and1111.

Ideally, in a perfect membrane state, stresses and strains would be distributed homogeneously across the thickness, and so would the growth multiplier ϑg. For this ideal state, we could have modeled skin as a perfect membrane. While the center of the expanded area is close to this ideal membrane state, the boundary experiences shear, and thus inhomogeneities across the thickness. Numerical investigations confirm that this boundary layer is restricted to approximately three rows of elements. We believe this is realistic and in good agreement with the boundary layer around the subcutaneous pocket, for example illustrated around the cheek expander in Figure 1.

Figure 12 shows that under the same pressure applied to the same base surface area, the circular expander in rows 1 and 2 induces the largest amount of growth, while the rectangular and crescent-shaped expanders placed orthogonal to Langer’s lines in rows 6 and 8 induce the smallest amount of growth. The contour plots confirm that it is important how the expander is positioned with respect to Langer’s lines. In summary, although the growth model itself is isotropic in the skin plane, the overall skin gain is larger when the expanders are placed with their long axis along the material’s strong direction, i.e., in the direction of Langer’s lines.

In this manuscript, we present a novel computational model for skin that allows to predict its acute and chronic behavior in respsonse to mechanical stretch within and beyond its physiological limit. The irreversible nature of the model allows to predict the long-term outcome of tissue expansion. Tissue expansion is a common procedure in reconstructive surgery that enables the body to grow extra skin to resurface large congenital defects, to reconstruct cancerous breasts, and to correct burn injuries. In tissue expansion, an inflatable silicone expander is implanted underneath the skin and gradually filled with saline solution. When skin is stretched beyond its physiological limits, new cells form and the skin grows. Despite intense research in skin growth, our understanding of the mechanobiological phenomena during tissue expansion remains poor and largely qualitative.

In this manuscript, we propose to model tissue expansion using the concept of finite growth based on the multiplicative decomposition of the deformation gradient into an elastic and a growth part. We assume that growth is an irreversible, transversely isotropic process which takes place exclusively in the skin layer, while the skin thickness is assumed to remain virtually unchanged. To model the chronic increase of skin area, we introduce a single scalar variable, the growth multiplier. Following clinical observations, we suggest that changes in this variable are governed by a stretch driven growth law. However, growth is activated only when skin is stretched beyond its physiological limit. Within the physiological range, we characterize the reversible behavior of skin through a transversely isotropic chain network model in which eight representative chains are modeled as a wormlike chains. The resulting elastic response captures the characteristic features of skin including its extreme nonlinearity, its locking behavior, and its anisotropic nature manifesting itself macroscopically through Langer’s lines. Since skin growth is a highly nonlinear, heterogeneous process, we propose to solve the governing equations using a nonlinear finite element approach. To model the chronic increase in skin area, we introduce a scalar-valued growth multiplier which we treat as an internal variable on the integration point level. We evaluate its temporal evolution locally using a finite difference approach. To guarantee maximum efficiency, stability, and optimal convergence of the algorithm, we suggest a local Newton iteration to update the growth multiplier at each integration point. Once we know the current amount of area growth, we evaluate the remaining elastic response using the wormlike chain network model. To update the deformation at each finite element node, we propose a global Newton iteration. This approach requires a consistent linearization of the biological equilibrium equation on the integration point level embedded within a consistent linearization of the mechanical equilibrium equation on the node point level.

To explore the conceptual features of our growth model, we study two simple model problems, displacement- and force-driven skin expansion. Upon displacement control, we observe material relaxation indicated through a gradual increase of area growth at a constant total stretch, while the elastic stretch and the resulting stresses decrease. Upon force control, we observe material creep indicated through a gradual increase of growth and of the total stretch, while the elastic stretch and the stresses remain constant. Lastly, we explore different expander geometries and study the effect of anisotropy by comparing different expander placements with respect to Langers lines. Despite the postulated isotropic in-plane area growth, we observe a significant sensitivity with respect to the expander orientation. Skin gain is found to be substantially larger when the expanders are placed along the direction of Langers lines.

Beyond the fractional area gain and the expander volume, which are used as characteristic metrics in this manuscript, Figures 10 and ​and11,11, our model is also capable to provides information about the heterogeneity of the tissue expansion process, Figure 12. While this might be less relevant in the context of pediatric forehead reconstruction, it is a critical issue in breast reconstruction, where the grown tissue needs to adapt a certain form and shape [16]. Our clinical collaborators have expressed interest in our computational model to optimize the expansion procedure in the context of heterogeneous growth to create different skin forms and shapes [70]. However, of course, next to the plain mechanical factors, there are always biochemical factors that might influence the choice of expander type, filling volume, and timing.

Although the proposed model for skin growth represents a tremendous advancement over the generic isotropic model we have recently proposed [15], we would like to point out that some limitations remain. First, motivated by experimental observations, which report normal cell differentiation upon tissue expansion [71], we have assumed that the material microstructure v remains unaffected by the growth process, i.e., Fg=ϑgI+[1−ϑg]n0⊗n0. It would be relatively straightforward to model the growth process itself as anisotropic [23]. This could imply growth ϑ|| exclusively along Langer’s lines, Fg=I+[ϑ||−1]ν0||⊗ν0||, or major in-plane growth ϑ|| along Langer’s lines combined with minor in-plane growth ϑ⊥ orthogonal to Langer’s lines, Fg=ϑ||ν0||⊗ν0||+ϑ⊥ν0⊥⊗ν0⊥+n0⊗n0. Similarly, we could even introduce a progressive reorientation of the collagen network to allow the material to align with the maximum principal strains [28, 36, 45]. To truly account for the microstructural origin of growth, which is beyond the scope of this manuscript, we could elaborate the use of mixture theories, which would allow us to model the turnover of the individual constituents such as collagen and elastin throughout the growth process [1, 29].

Second, for the sake of simplicity, we have modeled the tissue expander only implicitly through controlling the applied pressure. In real tissue expansion, the external control parameter is the expander volume [41]. In line with our discussion in Section 4.1, this implies that our virtual tissue expansion displays creep under constant loading, while clinical tissue expansion might rather display relaxation under constant deformation.

Third, here, we have assumed that the expander is connected tightly to the expanded tissue, neglecting effects of interface sliding and shear [62]. However, this seems to be a reasonable first assumption, since most current expanders have well-designed textures to promote mild tissue in-growth, primarily to prevent expander migration [9]. To address these potential limitations, we are currently refining the elastic model, the growth model, and the boundary conditions, to render our future simulations more realistic, and place it on an idealized face [43].

Last, while our computational model seems well suited to provide qualitative guidelines and trends, at its present state, it is not recommended for quantitative statements. We believe that using the equations on nonlinear continuum mechanics represents a significant advancement over the current gold standard to predict tissue growth exclusively in terms of kinematic quantities [59, 69]. However, acute and chronic in vitro and in vivo experiments will need to be performed to truly calibrate the underlying material parameters, to potentially refine, and eventually fully validate the model.

In summary, we have presented a novel computational model to simulate the acute and chronic response of skin when stretched within its physiological limits and beyond. This model captures the basic features of skin including its extreme nonlinearity, its locking behavior, its anisotropic nature, and its ability to grow when exposed to chronically elevated stretch levels. A comprehensive understanding of the gradually evolving stress and strain fields in growing skin may help the surgeon to optimize clinical process parameters such as expander geometry, expander size, expander placement, and inflation timing. Ultimately, through inverse modeling, computational tools like ours have the potential to rationalize these parameters to obtain skin flaps of desired size and shape. Overall, we believe that predictive computational modeling might open new avenues in reconstructive surgery and enhance treatment for patients with birth defects, burn injuries, or breast tumor removal.

This material was supported by the Claudio X. Gonzalez Fellowship and by the Mexican National Council of Science and Technology Scholarship CVU 358668 awarded to Adrin Buganza Tepole, and by the National Science Foundation CAREER award CMMI-0952021 and the National Institutes of Health Grant U54 GM072970 to Ellen Kuhl.

Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.